% Script to construct high resolution spectra of geomagnetic data
% Highly deterministic signals in the period range from 2 minutes to 20
% minutes were found in the HON 1 data's spectrum. The questions asked. Are
% they processing artifact ? Local communication noise ? Natural ?

% Load data 
%clear
time_array_min = (datenum(1995,1,1,0,0:9468000-1,30)); % 1995 - 2012
%time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2010,12,31,23,59,30))';

a = length(time_array_min);
%ncfname =  '/Users/manojnair/data/obs_mag_data/HON_1995_2012.nc';
%ncfname =  '/Users/manojnair/data/obs_mag_data/NGK_1995_2012.nc';
%ncfname = '/Users/manojnair/data/obs_mag_data/Canada/Absolute_netCDF/CBB_1995_2012_Absolute.nc';
%ncfname = '/Users/manojnair/data/obs_mag_data/image/Absolute_netCDF/HOR_1995_2012_Absolute.nc';
%ncfname = '/Users/manojnair/data/obs_mag_data/GUA_1995_2012.nc';
%ncfname = '/Users/manojnair/data/obs_mag_data/ASC_1995_2012.nc';
ncfname = '/Users/manojnair/data/obs_mag_data/CBI_1995_2012.nc';
[x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(ncfname, 0, a, 1);


%remove data gaps by interpolation

L = isnan(y_data);
y_data_interp = interp1(time_array_min(~L), y_data(~L), time_array_min,'cubic',NaN);
L = isnan(y_data_interp);
y_data_interp(L) = [];

L = isnan(x_data);
x_data_interp = interp1(time_array_min(~L), x_data(~L), time_array_min,'cubic',NaN);
L = isnan(x_data_interp);
x_data_interp(L) = [];

L = isnan(z_data);
z_data_interp = interp1(time_array_min(~L), z_data(~L), time_array_min,'cubic',NaN);
L = isnan(z_data_interp);
z_data_interp(L) = [];


% Now calculate the power spectra
% 60 days of nfft windows becuase of a suggestion
% by DJ Thomson


[Pxx1,F1] = pwelch(detrend(x_data_interp), hanning(1440* 60 ), 1, 1440 * 60 ,1/60);
[Pxx2,F2] = pwelch(detrend(y_data_interp), hanning(1440* 60 ), 1, 1440 * 60 ,1/60);
[Pxx3,F3] = pwelch(detrend(z_data_interp(10 * 1440 * 60 : 15* 1440 * 60 )), hanning(1440* 60 ), 1, 1440 * 60 ,1/60);


%plot the data

subplot(311);
loglog(F1 * 1e6, abs(Pxx1),'b');
axis([-inf inf -inf inf])
set(gca,'FontSize',16);
legend('HON X');
xlabel('\mu Hz')
ylabel('PSD');
subplot(312);
loglog(F1 * 1e6, abs(Pxx2),'b');
axis([-inf inf -inf inf])
set(gca,'FontSize',16);
legend('HON y');
xlabel('\mu Hz')
ylabel('PSD');
subplot(313);
loglog(F1 * 1e6, abs(Pxx3),'b');
axis([-inf inf -inf inf])

set(gca,'FontSize',16);
legend('HON Z');
xlabel('\mu Hz')
ylabel('PSD');

% % plot a zoomed in version for Z

figure (2);

plot(F1 * 1e6, abs(Pxx3),'b');
axis([2000 inf 0.0001 100])
set(gca,'FontSize',16);
legend('HON Z');
xlabel('\mu Hz')
ylabel('PSD');


for i = 1:7,
text(dp(i).Position(1),dp(i).Position(2),sprintf('%4.0fmuHz',dp(i).Position(1)));
end;



%% load and plot the 1 seconds data fro Pacific Islands

data = load('/Users/manojnair/projects/tsunami/IPM_PPT_2010_0910_mean_removed/IPM.nomean');
[Pxx1_PPT,F2] = pwelch(detrend(data(1:length(data)/2,1)), hanning(3600 *24 ), 1, 3600 * 24 ,1);
[Pxx2_PPT,F2] = pwelch(detrend(data(length(data)/2:length(data),1)), hanning(3600 *24 ), 1, 3600 * 24 ,1);
loglog(1./F2, abs(Pxx1_PPT),'r');
hold on;
loglog(1./F2, abs(Pxx2_PPT),'b');
